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SUDDEN STRETCHING OF A FOUR-LAYERED COMPOSITE PLATE 


by 

G. C. Sih 

Institute of Fracture and Solid Mechanics 
Lehigh University 
Bethlehem, Pennsylvania 18015 

and 

E. P. Chen 
Sandia Laboratories 
Albuquerque, New Mexico 87115 

ABSTRACT 

A research effort primarily concerned with the understanding of laminated 
composite plates with cracks subjected to time-dependent extensional loads is 
reported here. When loads are applied suddenly to a laminate, waves are re- 
flected and refracted through the laminae and give rise to stresses and strains 
throughout the composite system. The process is three-dimensional in character 
and presents a formidable problem in the theory of elastodynamics, particularly 
in the presence of crack-like imperfections. 

An approximate theory of laminated plates is developed by assuming that the 
extensional and thickness mode of vibration are coupled. The mixed boundary 
value crack problem of a four-layered composite plate is solved. Dynamic stress 
intensity factors for a crack subjected to suddenly applied stress are found to 
vary as a function of time and depend on the material properties of the laminate. 
Stress intensification in the region near the crack front can be reduced by having 
the shear modulus of the inner layers to be larger than that of the outer layers. 


INTRODUCTION 


The current interest in laminates for structural application is associated 
with the high strength-to-weight ratio which can be developed in laminates. 

These laminates are generally composed of layers which have been reinforced by 
embedding unidirectional fibers. The layers are adhered to each other such that 
the fiber direction varies from one layer to the next in a previously determined 
manner. The freedom of choice for fiber orientation in the layers of the com- 
posite system enables the development of laminates with special preferential 
directional properties for particular applications. Because of this character- 
istic of fibrous composites, the employment of these systems rather than equiva- 
lent homogeneous members will be clearly advantageous in many applications. 

Because of the complicated internal structure of composite systems, stress 
analysis is much more difficult than for equivalent single-phase material. One 
fact which emerges very clearly from laminate studies is that the stress field 
in composite systems is truly three-dimensional in character. Thus, even the 
stress field in a symmetric laminate subjected to in-plane loading cannot be ac- 
curately modeled by standard two-dimensional methods of analysis. The previous 
work in this area further indicates that relatively little effort has been made 
to formulate laminate plate theories that can effectively solve for the redistribu- 
tion of stresses and strains due to the presence of mechanical imperfections such 
as cracks. 

One possible means of simplifying the three-dimensional equations of elas- 
ticity is to invoke the concept adopted in the formulation of plate theory. Ap- 
proximate stress and strain dependence on the plate thickness coordinate are as- 
sumed such that the governing differential equations possess only two independent 


space variables. In addition, special attention must be given to the state of 
affairs near the crack when formulating plate theories for analyzing crack prob- 
lems. With this in mind, Hartranft and Si h [1] developed an approximate three- 
dimensional theory for a single material plate containing a through crack. The 
condition of plane strain was preserved ahead of the crack as suggested by Sih 

[2] . This theory was later extended to laminates by Badaliance, Sih and Chen 

[3] to solve the problem of a through crack in a laminar plate subjected to in- 
plane loading. The through crack configuration represents a preliminary effort 
to model the damage of composite plates. Additional complications arise when 
the load is time dependent. These considerations will be taken into account in 
the development of a new dynamic theory of laminated composite plates subjected 
to extensional loads. 

This work is concerned with the formulation of a dynamic theory of laminated 
plates and reduces to that of Kane and Mind! in [4] for the single material plate. 
The idealized condition of stress and displacement continuity across the inter- 
face is replaced by assigning certain conditions of material nonhomogeneity in 
the thickness direction of the laminated plate as if it were a single layered 
nonhomogeneous plate. The nonhomogeneity is made equivalent to a symmetric lami- 
nate balanced with reference to its mid-plane. A through crack is assumed to 
exist in a four-layered laminate. Dynamic stress intensity factors are obtained 
for the case of a suddenly applied uniform in-plane loading and shown to vary 
as a function of time. Discussed are also the influence of material properties 
of the layers on the local stresses. 


BASIC FORMULATION 


The elastodynamic equations of generalized plane stress are adequate only 
if the frequency of vibration is lower than that of the first thickness mode and 
the wave length is large in comparison with the plate thickness. In other words, 
the coupling between extensional and thickness mode ef vibration can be ne- 
glected. When laminated composite plates are stressed dynamically, loads are 
transmitted through the laminae by the reflection of thickness refraction of 
stress waves. The mode of vibration cannot be ignored, particularly in the vi- 
cinity of a crack-like imperfection where the stress state acquires a three-di- 
mensional character. 

A dynamic laminate plate theory will be developed to solve the problem of a 
four-layered composite plate with a through crack subjected to a suddenly applied 
uniform extensional load. The theory is a generalization to that given by Kane 
and Mindlin [4] for a single layer plate in which the extensional and thickness 
mode rf vibration are assumed to be coupled. Accounted for is the lowest thick- 
ness-stretch mode such that the displacement is normal to the plate surface. 
Mindlin and Medick [5] have also considered a formulation in which the thickness- 
shear mode of vibration with displacement parallel to the plate surface is also 
included. The mid-plane of the plate is taken as the nodal plane of vibration. 

Consider a four-layered composite plate of thickness h as shown in Figure 1 
where each layer has the same thickness h/4. The two outer layers have material 
properties (^>^2) or ( ^2 w ^ e the two inner layers have material properties 
(y-pv-j) or ( A-j ,y.| ) . The Lame coefficients are denoted by Aj and y. (j = 1,2). 

The layers are stacked such that symmetry prevails across the mid-plane of the 
laminate composite. The crack of width 2a cuts through the entire thickness of 
the laminate. 
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The time-dependent displacement field is assumed to be 


u x = v x (x,y,t) 

v y = v y (x,y,t) (1) 

w z = -jp v z (x,y,t) 

It follows that the strain components can be written as 
av x 

e x (x.y,t) = ^ 

9v v 

e y (*,y,t) - 
£ z (x,y,t) = f v z 

3V 3V ^ 

Y xy (x,y,t) = ay" + iT 

2 Z 9V 

Y yz (x,y,t) = h"F" 

2z 9v z 

y zx (x,y,t) 

in which the transverse normal and shear strains are assumed to be linear in the 
thickness coordinate z. If each layer of the laminated composite plate is iso- 
tropic, then the following stress-strain relationships may be used: 

a x = (X+2p)e x + A(ey+K£ z ) 

°y = ( A+2y ) £ y + *(« x + ™ z > 

a 2 = (X+2p)ic 2 £ z + XK(e x +e y ) 
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T yz = PY yz 
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T zx = yY zx 


T xy " yY xy 


(3) 


The constant 

k = it/vT2 


(4) 


accounts for the coupling between the extensional and thickness mode of vibration. 
It is determined from the three-dimensional equations of elasticity. As in the 
development of plate theories, the resultant strain quantities (A x ) > (Ay) > 


...» (O (j = 1.2) will be defined: 
xy 


h/4 


C(A V ) , (A v ) , (A,)_. (A xy ) ] = F m [ m Cs x > £ v » e z »Y xv ]dz 


1 " 1 


[(A*.) » (A v ) » (A 7 )_, (A xv ) ] - ^ ^ U f /A [ e x’ e y’ e z’ Y xy^ z 


' 2 y 2 2 xy 2 


h/4 

h/2 

h^4 


'x’ y’ z’ 'xy- 


-h/4 


+ L [e x » e y> e z »Y X y^ dz> 


■h/2 


(5) 


[(A xz ) » ( A yz \ ^ " h* j,A ^ Y xz’ Y yz^ zdz 


h/4 


xz'i’ ' yz 


-h/4 


.y Z - 


C(0 . ( a vt)J = Tib" . L [ Y x 2 ,Y v 2 ^ zdz + [ Y xz’ Y vzJ zdz 


-h/4 


h/2 


xz 'o’ ' yz'« J 7h 


■h/2 


xz’ ’yz- 


h/4 


’xz’ 'yz- 


Substituting equations (2) into (5), it is found that 


(Aj = (A v ) = 


3v. 


x y 2 3X 


- 7 - 




ivt — 

** 1 


3V. 


= ~ ay 


{h z\ s (A Z } 2 = F v z 


9v x + !^ 


^ A xy^ = ^ A xy^ ay + ax 


( 6 ) 


2 3v z 

^ A xz^, = ^ A xz^ ’ W ax” 


1 


. 2 3v z 


(A yz\ = (A yz } 0 = F ay 


The laminate plate theory can be most conveniently formulated in terms of the 
stress resultants 


h/2 

(N„,N„,N,,N V „) • [_ (o x ,a y .o z ,t xy )dz 


x* y z xy 


-h/2 


(7) 


and the transverse shears 


h/2 

(R..R„) * L <V T yz )2dz 


x’ y 


■h/2 


(8) 


The stress-strain relations in equations (3) when enforced yield 


h av 9v 

N v (x,y,t) = j [(B+2y) -%r- + 3 gT^l + 3 kv. 


h av v av x 

N w (x,y,t) = j [(S+2y) + B -tH + 3kV 


3X J " ”Z 
9V.. 9V, 


2 • 9y 

N z (x,y,t) = (s+2 y )k 2 v z + ^ B<h (g~ + 

1 9v x 9v v 
N xy (x,y,t) = 2 Y h (ay- + aT 5 


( 9 ) 
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and 


R x (x,y,t) - lg sh 2 

1 8V 7 

R y (x»y,t) = 38 5h2 W 


(10) 


In equations (9) and (10), g, y and 6 stand for 


3 = x ] +;v 2» Y = V‘-| + v*2» 6 = P 1 + ^ y 2 (11) 

Denoting and p£ as the mass density of the inner and outer layers of the lami 

nate, the three equations of motion governing N , N R become 

x y y 


( 12 ) 

a5T T W ' "z " 48 


9N x 9N xv 1 92v x 

+ ~dy L = 2 h ( p l +p 2^ at 7- 


ax 


9 2 v. 


!!k + = i h(p +p , Lii 

3x + 3y 2 n ' p l p 2' 3t 2 


3R. 


9R y 


_ 1 


u2 ( _ _l“7 ^ \ 


The result of substituting equations (9) and (10) into (12) renders 


~ 8V 9V\. r\ Q a V 

yv 2 v x + (3+y) (p l +p 2 ) 


3V_ 


3 2 V, 


at 9 


. 3V 3V 3V 3 2 V 

W 2 v y + (6 +y) 3y (-aT- + if' 1 + TT 3T = (,, 1 +P 2 ) 3t^ 


(13) 


/IQ 9/1d 8V v 9V w 3 2 V, 

^V 2 v 2 - w (3+2 Y )k 2 v z - (gjf + -§y^) (p-,+7p 2 ) gp- 
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where v 2 = 3 2 /3x 2 + 3 2 /3y 2 is the Laplacian operator in two dimensions. 

METHOD OF SOLUTION 

Equations (13) will be solved by introducing two potential functions 4(x,y,t) 
and H(x,y,t) as 

v x( x,y ’t) = H + 9y 

(14) 

v y (x,y,t) ay " ax 

Making the appropriate algebraic manipulations, the governing equations for the 
potential functions can be derived by enforcing equations (13): 

yV 2 H = (p,+p 2 ) 0 - 6(8+2y)v"* + 192 

- <(p-|+P 2 )C(p-| + 7p2) ft^ + (S+Z'l) 

- [4(p,+p 2 ) + (B+2 yMp 1 +7p 2 ) 0 (v 2 *)]} (15) 

Once 4(x,y,t) and H(x,y,t) are known, v x and v y can be obtained from equations 
(14) and 

v z (x,y,t) = ^ C(p-i + p 2 ) ft^ " (S +2k ) v ^^ 

Suppose that a uniform stress resultant N Q is applied suddenly to the crack 
surfaces and the resulting deformation is symmetrical with respect to the x-axis 

then the following' conditions are to be specified: 
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N y (x,o,t) = -N Q H(t), x<a 

07 ) 

v y (x,o,t) = 0, x>a 

where H(t) is the Heaviside unit step function. The condition of symmetry fur- 
ther requires that 


N xy (x,o,t) = R y (x,o,t) = 0, for all x (18) 

'At * 

Use will now be made of the Laplace transform. Let <p (x,y, p), H (x,y,p), etc., 
denote the Laplace transforms of the functions <f>(x,y,t) , H(x,y,t), etc. Equa- 
tion (15) when expressed in the Laplace transform domain become 

(v 2 -o)|)fJ(x,y,p) = 0 

(v 2 -u|)<|)2(x,y,p) = 0 (19) 

(v 2 -to|)H*(x,y,p) = 0 

where the potential <}>(x,y,t) has been separated into two parts: 

<i>(x,y,t) = 4 >-|(x,y,t) + <j> 2 (x,y,t) (20) 

in terms of time t or 

it it it 

<f» (x,y,p) = <j» 1 (x,y,p) + $ 2 (x,y,p) (21) 

in terms of the Laplace transform variable p. The parameters coj (j = 1,2,3) in 
equations (19) are defined as 
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( 22 ) 



<•»§ ■ (p-|+P 2 )P 2 y " 1 


in which the newly defined quantities are 

“o ” (p-j+p^Yq’ 5 o (p 1 +7p 2 )yq 

_ 2 - 12 k 2 (8+2y) 

p o Pi+7p 2 * o h ii (p^+p 2 ) 

and y q takes the form 

y 2 = 4y(e+y) 

Y 0 (p-|+p 2 )(e+2y) 

Equations (19) then give 


fj(x,y,p) = 7 f A(s,p) cos(sx) exp(-s 1 y)ds 


4>o(x,y,p) =f / B(s,p) cos(sx) e?<p(-s 2 y)ds 

<i ir o 


H*(x,y,p) = 7 / C(s,p) sin(sx) exp(-s^y)ds 
* o ■* 


with s. being given by 

vJ 


s. = /i^, j = 1,2,3 


(23) 


(24) 


(25) 


(26) 
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The dynamic problem has now been reduced to finding the three unknown functions 
A(s,p) , B(s,p) and C(s,p). 

DUAL COUPLED INTEGRAL EQUATIONS 


Before the boundary and symmetry conditions can be enforced, it is necessary 

He He 

to obtain v ¥ (s,y,p), v (x,y,p), etc., in terms of the unknowns in equations 
x y 

(25). With the help of equations (14) and (16), it can be shown that 


v v (x,y,p) = - 7 / [sA(s,p) exp(-s,y) + sB(s,p) exp(-s«y) 
* o 

+ s 3 C(s,p) exp(-s 3 y)] sin(sx)ds 


v*(x,y,p) = - 7 / [s.,A(s,p) exp(-s,y) + s 2 B(s,p) exp(-s 2 y) 

y O' 


(27) 


+ sC(s,p) exp(-s 3 y)] cos(sx)ds 


v*(x,y,p) =- / Ca 1 A(s,p) exp(-s 1 y) + A 2 B(s,p) exp(-s 2 y)] cos(sx)ds 


The quantities a. (j = 1,2) are given by 

J 

. . h(e+ 2 y) r (p i +p 2 )p2 21 i - 1 ? 

A j " 26 k 0+2 y " J 1,2 


(28) 


He He 

Similarly, the Laplace transform of N x (x,y,p), N y (x,y,p) become 


* 9 (p,+P9)p 2 

N x (x,y,p) = 7 yh / ([ ^ s 2 ] A(s,p) exp(-s 1 y) 


(Pi + Po)P 2 

+ [ — s|] B(s,p) exp(-s 2 y) - ss 3 C(s,p) exp(-s 3 y) }cos(sx)ds 
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* 2 °° (Pi + P?)p 2 

N y (x,y,p) = - yh / {[s 2 + ^ ][A(s,p) exp(-s.jy) + B(s,p) exp(-s 2 y)] 

+ ss 3 C(s,p) exp(-s 3 y)}cos(sx)ds 

(29) 

N xy (x,y,p) = - yh / {ss-|A(s,p) exp(-s-|y) + ss 2 B(s,p) exp(-s 2 y) 


+ j (s 2 +s|) C(s,p) exp(-s 3 y)}sin(sx)ds 
while R x (x,y,p) and R y (x,y,p) take the forms 

* jjL2 00 

R x (x,y,p) = - 2^ / s[A-|A(s,p) exp(-s,y) + A 2 B(s,p) exp(-s ? y)] sin(sx)ds 

0 t c c 

m 

Ry(x,y,p) = - Ifc" / Cs-|A-jA(s,p) expC-s^) * s 2 A 2 B(s,p) exp(-s 2 y)] cos(sx)ds 


The symmetry conditions in equations (18) when applied show that A(s,p), B(s,p) 
and C(s,p) can be expressed in terms of a single unknown D(s,p): 


A(s,p) = 


S 2 + S 2 


D(s,p) 


B(s,p) = -s 1 


r (p-j+P 2 )p 2 

^ 6+2y 


0)|]/{S 2 


[■ 


(p1+P 2 )p 2 

3+2y 


- u|]}A(s,p) 


(31) 


C(s,p) = 


2ss n (to|- 


■»|) 


(Pi + Po)P Z 

' s2+s §)c-y^ 


u,3 


A(s,p) 


Application of the mixed boundary conditions in equations (17) leads to a sys- 
tem of dual integral equations 


/ D(s,p) cos(sx)ds = 0, x>a 
o 

» irN 

/ sF(s,p) D(s,p) cos(sx)ds = - x<a 

0 T K 


( 32 ) 


The function F(s,p) is known: 


F(s,p) = 


s2 + s § 

ss. 


CCS 2 + 


(p-| + P2)P 2 

2y 


•] {1 - 


S 1 [ 


(Pi +P2 ) P 2 


J+2y 


- to 2 ] 


2s2s 1 s 3 


2 2 
o)|-o)2 


s 2 [ 


(p-i+p 2 )p 2 

3+2y 


to 


§] 


[ ( Pl +p 2 ) p^/ ( e+2r 


(33) 


The standard procedure by Copson [6] may be applied to solve equations (32) and 
the result is 


D(s,p) = 


"N a 2 C(p-| + Pp)P 2 /( 3 + 2y) ]- oj| 

W 


-} 


(1 - 


e+57)( (0 i' a> 2) 


x / /£ $*U,p) J (sac)ds 

0 0 


(34) 


★ 

in which $ U,p) can be computed from a Fredholm integral equation of the second 
kind: 


ic ■ * 

$ U>P) + f <t> (n.p) L(s»n>p)dn = 
0 

The kernel L(€,n>p) is 


(35) 


L(c>n>p) *^7 s[G(|, p) - 1] J 0 (sg) J Q (sn)ds (36) 

o 
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while the function G(s,p) Is related to F(s,p) in equation (33): 




{[(p-i + P2)P 2 /(6+2y)]-oj|}y 
G(S ’ P) = P z Cl-Y/(s+2 Y )Ku,f-oi|) F(s - P> 


DYNAMIC STRESS INTENSITY FACTORS 


( 37 ) 


Of interest is the intensification of the dynamic stresses ahead of the crack. 
Hence, the integrals in equations (29) and (30) must be evaluated for large values 
of s which corresponds to distances near the crack edge x = ±a and y=0. In terms 
of the polar coordinates r and 6 in Figure 1, the Laplace transform of the stress 
resultants for small r are found: 


N*(r,e,p) = jjfi cos i (i _ sin | sin |^-) + 0(r°) 

/2r 

k*( 

N*(r,e,p) = - — cos 4 (1 + sin 4 sin P-) + 0(r°) 

y v^r d d 2 

(38) 

N* (r, Q ,p) = — sin 4 cos 4 cos P- + 0(r°) 

xy c c c 


R*(r,e,p) = R*(r,e,p) = 0(r°) 
in which k*(p) is the Laplace transform of k-j(t): 
k*(p) • N 0 /f 


( 39 ) 


The Laplace inversion theorem may now be applied to give 
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k (t) 

N v (r,9,t) = — cos I- (1 - sin sin |^) + 0(r°) 

x v^F 2 2 2 

N w (r,0,t) = — cos (1 + sin f- sin -P-) + 0(r°) 

y /ft 2 2 2 

Mt) fi fl n 

N v ,.(r,0,t) = sin ■§• cos I- cos 4r- + 0(r°) 

xy c c c 

R x (r,0,t) = R y (r,0,t) = 0(r°) 


( 40 ) 


Equations (40) reveal that dynamic loading does not affect the functional re- 
lationship of r and 0. The stress intensity factor, however, is a function of 
time: 

N.t/a * n * 

^(t)*^-/ — exp(pt)dp (41) 

★ 

where Br denotes the Bromwich path of integration. Once $ (5,p) is calculated 
from equation (35) and evaluated at £= 1, equation (41) may be solved numerically. 

* vl/2 

Figure 2 gives a plot of $ (l,p) as a function of C2-|/pa where C21 = (u-j/p-j) 
is the shear wave velocity referred to the material in the inner layers. For 

•k 

P 1 = p 2’ V 1 = v 2 = ^*3 anc * P 1 = P 2 J $ 0»P) is seen t0 increase monotonically 
with C 2 -j/pa* Three different ratios of ^ 2/^1 = 0*2, 1.0 and 5.0 are considered. 
Making use of the results in Figure 2, k^(t) in equation (41) may be computed. 

Refer to Figure 3 for a display of k.|(t)/N 0 vef versus C 2 -jt/a. The resultant 
stress intensity factors are observed to vary as a function of time. Their am- 
plitude rise quickly reaching a peak and then declines. The solution for a homo- 
geneous plate corresponds to y 2 /y-| = 1.0 as the Poisson's ratio and mass density 
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0 12 3 4 5 


pa 

Figure 2 - Numerical values of $ (1 , p) as a function of c ?1 /pa for a/h 



for the inner and outer layers are assumed to be equal. The peak value of k-j ( t ) 
is greater than that of the homogeneous plate solution for yg^-j while the opposite 
is found for yg^-j. Hence, the intensity of the crack border stress field can 
be reduced by having the shear modulus of the outer layers to be smaller than 
that of the inner layers. 

CONCLUDING REMARKS 

A dynamic laminate plate theory has been developed for solving crack bound- 
ary value problems. The complexity of the problem owing to material nonhomogeneity 
and dynamic stress analysis necessitates certain simplifying assumptions so that 
effective analytical solutions can be obtained. It is shown that the dynamic 
stresses near a mechanical imperfection such as a crack are intensified depending 
on the stacking sequence of the laminae. In general, this intensity tends to in- 
crease quickly for small time reaching a peak and then decreases to the static 
solution for sufficiently long time. When the modulus of the outer layers are 
smaller than that of the inner layers, the crack border stress intensity reaches 
a maximum quicker than the homogeneous solution but with a smaller magnitude. 

The opposite holds for the case when the outer layers are stiffer than the in- 
ner layers. Information of this type is useful for evaluating the resistance of 
laminate plates to impact loadings. 
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COMPUTER PROGRAM: DYNAMIC LAMINATE PLATE THEORY WITH A CRACK 


1 


PROGRAM FLAP (INPUT. OUTPUT) 

REAL NON (4) *F(4.4»?).G(4*4>*D<4).PT 
REAL 0(4) *C (4 ) 

REAL LP (19) *OTA( 19) 

5 


EQUIVALENCE (NONtR) 

COMMON KltK2*K3.K4 
COMMON/AUX/H»P»PKl ,PK2*BMU» X»Y 
LP(1>=0.0 
OTA ( l ) =0 • 0 

10 


READ 2*K1 *K2*K3»K4 
2 FORMAT (12) 



K 1 = ORDER OF SYSTEM OF EQUATIONS 


• 

K2 * NO. OF DISTINCT KERNELS 


* 

K3 = NO. OF DATA POINTS 

15 


K4 * NO. OF DATA SETS TO BE EVALUATED 


* 

SET UP DATA POINTS 
AK=K3 

DO 5 N=1 *K3 

an=n 

20 


5 PT(N)=AN/AK 


• 

SET UP INTEGRATION MATRIX 
M=K3-2 
N=K3"1 
A =K 3 

25 


A = l ,/(3.*A) 

DO 10 K=2 *M « 2 
10 0(K)=2.*A 

DO 15 K = 1 *N » 2 
15 0 ( K ) =4 • * A 

30 


0(K3)=A 


* 

CALCULATE NOMHOMOGENFoUS TERMS 
PHS=1 .0 
DO ?2 1=1* K2 
PRINT 9 

35 


9 FORMAT ( 1 H 1 ) 

DO 999 11 = 1 » K4 
DO 35 N= 1 » K 3 

35 NON(N)=RHS*SQRT (PT (N) ) 


CALL CONST (I) 

* CALCULATE kernel matptcfs 
00 20 N=1 *K3 
DO 20 M = 1 ♦ K 3 

F ( M *N * I ) =EU ( I .PT (M) ,PT (N) ) 
20 CONTINUE 

CALL CHANGE(F»G.D. T ) 

CALL LINEO (G*B*C* K 3 ) 

DO 40 L=1*K3 
PRINT G*PT <L) *NON U. ) 

6 FORMAT (5X*F8»4«F15.6) 

4 0 CONTi 'UE 

L P < 1 1 ♦ 1 ) =NON ( K 3 ) 
DTA(II*1)=p 
999 CONTINUE 

CALL LAPINV(OTA.LP) 

2? CONTINUE 
FNO 


45 


50 


55 
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1 


5 


10 


IS 


20 


?S 


30 


35 


40 


FUNCTION SIMP ( I *A«P) 

C0MM0N/AUX/H»P»PK1 «PK2 ♦ 8MU ♦ X » Y 

MXVZs2**15 

DEL»0 *25* (B-A) 

IF (DEL ) 40 »45 *50 
45 SIMP=0*0 
RETURN 
50 CONTINUE 

SArZUtA) *Z ( I ♦B) 

SB*Z ( I * A *2 * •DEL ) 

SC=Z ( I ♦A*OEL) ♦ Z < I ♦ A*3«*0EL> 

Si* (OEL/3 • )*(SA»2»*SR^4**SC> 

IF(SI.EO.O.O) GO TO 45 
K = 8 

35 SB=SB*SC 

DEL=0.5*OEL 
SC=Z U«A*DEL) 

J“K - 1 

00 5 N=3.J»2 
AN=N 

5 SC=SC*Z(I»A»AN»DEU 

S2= (DEL/3. )<MSA*2.*SR*4.*SC) 

OlF = ABS ( (S2-SD/S1) 

ER=0.01 

IF (OlF-ER) 30*25.25 
30 SIMP=S2 

return 

25 K = 2*K 
Sl = S2 

IF (K-MXYZ) 35.35.40 
40 PRINT 42 ♦ I « A » B 

42 F0RMAT(5X«* INT. DOES NOT CONVERGE *»I3*2F9.4) 
PRINT 60 « X * Y 
60 FORMAT (2F10. 5) 

00 70 J* 1 ♦ 1 0 
D I P = J 

OIP=OIP/10. 

W=Z(I»OIP) 

PRINT 60 1 W 
70 CONTINUE 
CALL EXIT 
END 


SYMBOLIC REFERENCE MAP <R=1> 

ENTRY POINTS 
4 SIMP 


REAL 

PEAL 

REAL 

REAL 


variables 

SN TYPE 

0 

A 

REAL 

0 

B 

REAL 

2S0 

del 

REAL 

264 

DIP 

REAL 


relocation 


F.P. 

260 

AN 

F.P. 

4 

Bmu 


2 62 

OIF 


263 

ER 
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SUBROUTINE CHANGE (c* »G*0» I ) 
COMMON K 1 1 K2 «K3 1 K4 



REAL E (A»4*2) « G ( 4 . 4 ) » D (4 ) 



DO 

10 N=1*K3 



5 

DO 

10 M=1*K3 




G<M 

,N) =F(M*N*I)*n<N) 




10 CONTINUE 




DO 

20 N=1«K3 




20 G(N 

«N) =G (N«N ) ♦ 1 #0 



10 

RETURN 



■ 

END 




SYMBOLIC 

REFERENCE 

MAP ( R= 1 ) 



ENTRY POINTS 





3 CHANGE 





VARIABLES sn 

l TYPE 

relocation 



0 0 

REAL 

ARRAY F.P. 

0 

E 

0 G 

REAL 

ARRAY F.P. 

0 

I 

0 K 1 

INTEGER 

/ / 

1 

K2 

2 K 3 

INTEGER 

/ / 

3 

K4 

53 M 

INTEGER 


52 

N 

STATEMENT LABELS 





0 10 


0 

20 


LOOPS LABEL 

INDEX 

FROM-TO LENGTH 

PROPERTIES 

17 10 * 

N 

4 7 17B 


NOT 

70 10 

M 

5 7 7R 

INSTACK 


43 20 

N 

B 9 4 B 

INSTACK 


COMMON blocks 

LENGTH 




/ / 

4 





REAL 

INTEGER 

INTEGER 

integer 

INTEGER 


JNNEP 


STATISTICS 

PROGRAM LENGTH <S5B 53 

SCM BLANK C.GMMON LENGTH 4B 4 

47000R SCM USED 
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if— ' TTET.j 


1 


5 


10 


15 


20 


25 


30 


35 


40 


SUBROUTINE LINEQ(A.B.T.N) 
PEAL A (N.N) «B(N) «T(N) 

DO S I =2 ♦ N 

5 A(Ifl)sA(I«l)/A(l«I> 

DO 10 K=2»N 
M=K«1 

DO 15 1=1 *N 
15 T ( I ) =A ( I «K) 

DO 20 J=1.M 
A < J * K ) =T ( J) 

Jl=JM 

DO 20 I=J1.N 
T<I)=T(I>-A(I.J)»A(J,K) 

20 CONTINUE 

A <K,K)=T (K) 

IF(K.EO.N) GO TO 10 
M=K* 1 

DO 25 I =M ♦ N 
25 A(I.K)sTfI)/A(K,K) 

10 CONTINUE 
* BACK SUBSTITUTE 
DO 31 1=1. N 
T ( I ) =B ( I ) 

M= I ♦ 1 

IF(M.GT.N) GO TO 31 
DO 30 J=M.N 
8(J)=B(J)-A(J,I)*T(I) 

30 CONTINUE 
31 CONTINUE 
DO 35 1=1. N 
K=N ♦ 1 - I 

P(K)=T (K)/A(K.K) 

K 1 =K - 1 

IF(KI.EO.O) GO TO ?5 
DO 36 J1 = 1 * K 1 
J=K-J1 

T(J)=T<J)-A( J.K)*R(K) 

36 CONTINUE 
35 CONTINUE 
PETUBN 
END 


SYMBOLIC REFERENCE MAP (R=l) 


ENTRY POINTS 
3 LINED 


variables 

0 A 

SN TYPE 
REAL 

RELOCATION 
ARRAY F.P. 

0 

B 

PEAL 

172 

I 

INTEGER 


1 75 

J 

INTEGER 

176 

J1 

INTEGER 


173 

K 

INTEGER 

177 

K 1 

INTEGER 


174 

M 

INTEGER 

0 

N 

INTEGER 

F.P. 

0 

T 

REAL 
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1 


5 


10 


15 


FUNCTION FU { I * A *B) 
COMMON/AUX/H.P*PKl ,PK2.RMU*X.Y 
X = A 
Y=B 

IF ( A*B ) 5*10*5 
10 FUsO.O 
RETURN 

5 SUM=SIMP<I»0. 0*5.0) 

ER*0.01 
DEL *5.0 
20 UP*DEL*5.0 

ADDL*S IMP (I » DEL .UP ) 

DEL =UP 

TEST=ABS (ADDL/SljM) 

SUM=SUM+ADDL 
IF(TEST-ER) 15.20*20 
15 FU=SORT (X»Y)«SUM 
RETURN 
END 



symbolic 

REFERENCE 

MAP 

( R=1 ) 




ENTPY 

POINTS 







4 

FU 







variarles SN 

TYPE 


RELOCATION 




0 

A 

REAL 


F.P. 

62 

ADDL 

real 

0 

B 

PEAL 


r.R. 

4 

BMU 

REAL 

60 

DEL 

REAL 



57 

FR 

REAL 

55 

FU 

REAL 



0 

H 

REAL 

0 

I 

INTEGER 


F.P. 

1 

P 

REAL 

2 

PK 1 

REAL 


AUX 

3 

PK2 

REAL 

56 

SUM 

REAL 



63 

test 

real 

61 

UP 

REAL 



5 

X 

real 

6 

Y 

REAL 


AUX 




externals 

TYPE 

APGS 






SIMP 

REAL 

3 



SORT 

real 

INLINE 

FUNCTIONS 

TYPE 

ARCS 






ABS 

REAL 

1 

INTRIN 





STATEMENT LABELS 
1 A 5 
2 ? 20 


10 


INACTIVE 


COMMON BLOCKS LENGTH 
AUX 7 


STATISTICS 

program length 

SCM LA8ELED COMMON LENGTH 
470008 SCM USED 


64B 

7R 


52 

7 
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1 FUNCTION BESJO ( A ) 

IE(A-3.)5*5*10 
5 B=A*A/9. 

Vf»l. -2. 2499997*8 
5 Z=B*B 

W*W4l.2656208*Z 

Z=Z*8 

W*W-*3163866*Z 

Z S Z*B 

10 W*W*.0444479*Z 

Z s 7*8 

W S W-* 00 39444* Z 
7=Z*B 

BESjOaW* *00021 *Z 
15 RETURN 

10 B=3./A 

W=*79788456-.OOOOOft77»B 
V=A-,78539816-*041 *,6397*8 
Z=B*8 

20 W S W- • 00552 74* Z 

V S V«* 000 03954* Z 
Z=Z*B 

W*W-.00009512*Z 
V*V^ *0026257 3*Z 
25 Z=7*B 

W*W ♦,001372 37 *Z 
V*V-*00054125*Z 
Z=7*B 

W»W-,00072805*Z 
30 V*V-. 00029333*2 

Z*Z»B 

W=VU. 00014476*7 
V=V>*00013558*2 
BESjO*W/SQRT (A)*CO« » V > 

35 RETURN 

END 


SYMBOLIC REFERENCE MAP <R=1> 

entpy points 

4 BESJO 


VART ARLES 

SN TYPE 

RELOCATION 




v 0 A 

REAL 

F.P. 

114 

8 

real 

113 BESJO 

PEAL 


117 

V 

real 

11B W 

PEAL 


116 

z 

REAL 

EXTERNALS 

TYPE 

ARGS 




COS 

REAL 

1 LIBRARY 


SORT 

real 


STATEMENT LABELS 
0 5 


INACTIVE 


26 10 


1 


5 


10 


SUBROUTINE CONST ( I ) 

COMMON/ AUX/H* P*PK1 .PK2»BMU*X ? Y 
PKlsO.3 
PK2=0.3 
BMU=50.0 
H=1 .0 
READ 2 *P 
2 FORMAT (F10. 5) 

HH=1 ./H 

PRINT 1»BMU.PK1,PK?,HH,p 

l FORMAT (/////5X** mU?/MU1 =*F6.2** NU1 s*F4.2«* NU2 « 
1 A/H =*F4.2.* C2 l/°4 =*F4.2//) 

RETURN 

END 


SYMBOLIC reference MAP <R=1) 


entry points 
3 CONST 


variables 

SN TYPE 

relocation 




4 

BMiJ 

real 

AUX 

0 

H 

real 

55 

HH 

REAL 


0 

I 

INTEGER 

1 

P 

REAL 

AUX 

2 

PK 1 

REAL 

3 

PK 2 

REAL 

AUX 

5 

X 

REAL 

6 

Y 

real 

AUX 





FILF NAMES MODE 

INPUT FMT OUTPUT FMT 


STATEMENT labels 

37 1 FMT 25 2 FMT 

COMMON BLOCKS LENGTH 
AUX 7 

STATISTICS 

PROGRAM LENGTH 560 46 

SCM LABELED COMMON LENGTH 7B 7 

470008 SCM USED 
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1 


5 


10 


15 


?0 


25 


30 


35 


40 


FUNCTION Z(I*S) 

COMMON/ AUX/H»P*PK1 ,PK2,9MU*X*Y 
COMPLEX DA*DL1*0L2.SA.SB*SC*SD 
COMPLEX GA* G8 * CA ♦ CR * CC ♦ F * G 
PI*3. 1415926 
PP=P*P 

PG*2*/PP/< 1 •♦BMW) 

AAs2.*0 # -PK1)/<1.-2.*PK1) 

ABs2**(l.-PK2)/(t.-2,*PK2) 

PA*2,/PP/ ( AA*9MU* AR ) 

PO*2#*H*H/PI/PI/ (AA*BMU*AB)/PP 
PA= ( 1 . ♦BMU) / ( AA*8Mi i*AB ) 

BB=1 .-BA 

BCa(1.^7.*BMU)/4./(l .♦BMU) 

B0*Pl*PI/2./H/H 
ALP*1 ./4a/BA/8B 
DLPs9C/4./PB 

DO*( <ALP*OLP)*POM ,)**?-4.*ALP*0LP*P0*(P04l 
G = CMPLX (00*0.0) 

OAsCSORT(G) 

DLlsBD/DLP* ( ( ALP^Ol.P ) *PO* 1 • ♦DA) 

DL2sBO/OLP* ( (ALP*0|_P) ttp O*l.-DA) 

SC=S*S*DL1 

SDsS*S*DL2 

GAsCSORT(SC) 

GB=CSQRT ( SO ) 

GCaSORT (S*S*PG) 

SAa(PA-0L2)/(0Ll-DL2) 

SBa (PA-DLl ) / (PA-DL 7 ) 

CAaSA/PG/BB 

CBa?.*(S*S*PG/2. ) *»2/GA* ( 1 .-GA/G8*SB) 

CC=2.*S*S*GC/SA 

F=CA*(C3-CC) 

Q=PEAL (F) 

OAaAlMAG(F) 

IF(OA-O.0)5» 10*5 

10 Z=(Q-S)»BESJO<S*X)*RPSJO<S*Y) 

PETURN 

5 PRINT 9*P*S*F 
9 FORMAT < 4F 10.5) 

CALL EXIT 

end 


SYMBOLIC REFERENCE MAP (Pal) 

entry points 

4 Z 


variables 

SN TYPE 

RELOCATION 



276 

AA 

REAL 


277 

AB 

306 

ALP 

REAL 


302 

BA 

303 

BB 

REAL 


304 

BC 

305 

BD 

REAL 


4 

BMU 


REAL 

real 

real 

PEAL 
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1 SUBROUTINE LAPINV <OL AM .PHI ) 

C THIS PROGRAM EVALUATES THE COEFFICIENTS FOR SERIES 

C OE JACOBI POLYNOMIALS WHICH REPRESENTS A LAPLACE 

C INVERSION INTEGRAL 

5 REAL MUL 

DIMENSION A (50) ♦GLAM (50) .PHI (50) .C<4.50> 

DIMENSION 9K(101I .TT(101) 

C0MM0N/2/T1 *TF»DT.MN»8K.TT 
READ l.NN.MN.MM 
10 1 FORMAT (312) 

READ 2*TI *TF »DT 

2 FORMAT (3F10.5) 

PRINT 99 

99 FORMAT < 1H1 ) 

15 CALL SPLICE (GLAM.PhI.MM.C) 

PRINT 101 

101 FORMAT (/////5X.* GLAM PHI *) 

PRINT 102* (GLAM ( 1 ) ,PHl ( I ) * 1*1 *MM) 

102 FORMAT (5X*F10*5.5X.F10*5) 

20 Ml 1=MM-1 

PRINT 99 
DO 10 1*1 *NN 
RE AO 3.BET .DEL 

3 FORMAT (2F10.5) 

25 PRINT 98.BET .DEL 

98 FORMAT (/////5X.*8FtA *»F5.3.» DELTA =*F5.3) 

DO 11 L=1*MN 
AL*L 

S=1./(AL»BET)/DEL 

30 CALL SPLINE (GLAM*P-iI .MM.C.S.G) 

F=G*S 

IF (AL-2.181 .82.83 
81 A(l) = <l.*BET)*DEL*c- 
GO TO 11 

35 82 A(2)*(<2.*BET)*0EL*F-AU))M3.*BET) 

GO TO 11 
83 CONTINUE 
TOP* 1 • 

L 1 =L” 1 

40 AL 1 =L 1 

00 12 J*1 *L1 
A J= J 

TOP=AJ«TOP 
12 CONTINUE 
45 L2=2*L“1 

BOT*l • 

DO 13 J=L.L2 
A J* J 

eOT=(AJ*PET)*BOT 
50 13 CONTINUE 

MUL=BOT/TOP 
SUM=0.0 
DO 14 N=1 *L1 
AN=N 

55 IFUN-2. 185*86. 87 

85 TOD=l. 

GO TO P8 


_4f 
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60 


06 TOO=ALl 
GO' TO 88 

87 CONTINUE 
TOD=l • 
lCH=Ll-(N-2) 

DO 15 J=ICH*L1 
AJ- J 

65 TOD=AJ*TOD 

15 CONTINUE 

88 CONTINUE 
BODsl. 
ja=li*n 

70 DO 16 J=LtJA 

A J= J 

BODsBOD* ( A J*8ET ) 

16 CONTINUE 
C0=T00/80D 

75 SUMsSUM*CO»A <N) 

14 CONTINUE 

A (L)=MUL*<DEI*F-SUm) 

11 CONTINUE 

CALL JACSER(DEL«A,cET> 
80 10 CONTINUE 

999 CONTINUE 
RETURN 

end 


symbolic reference 


ENToy POINTS 
3 LAPINV 


VARIABLES 

SN TYPE 

377 

A 

REAL 

354 

AL 

REAL 

371 

AN 

REAL 

4 

8K 

REAL 

366 

80T 

REAL 

376 

CO 

REAL 

2 

DT 

REAL 

366 

G 

REAL 

347 

I 

INTEGER 

363 

J 

INTEGER 

353 

L 

INTEGER 

365 

L2 

INTEGER 

3 

MN 

INTEGER 

350 

Mil 

INTEGER 

3^ 

NN 

INTEGER 

355 

S 

REAL 

1 

TF 

REAL 

37? 

TOO 

REAL 

1 c; i 

TT 

ftp AL 


MAP <R=1 ) 

RELOCATION 

ARRAY 

ARRAY 2 
2 

2 

2 


364 

AJ 

real 

362 

AL 1 

real 

351 

BET 

REAL 

374 

800 

real 

461 

C 

real 

352 

DEL 

real 

357 

F 

real 

0 

GLAM 

REAL 

373 

ICH 

INTEGER 

375 

JA 

INTEGER 

361 

LI 

integer 

346 

MM 

INTEGER 

344 

MUL 

REAL 

370 

N 

INTEGER 

0 

PHI 

real 

367 

SUM 

real 

0 

T I 

real 

360 

TOP 

real 


ARRAY 


2 




SUBROUTINE JACSER (n*C*B) 

DIMENSION C(50) fSF(50) fP(50) 

DIMENSION BK(lOl) .TT(lOl) 

C0MM0N/2/TIfTFfDTfMNf8KfTT 

TT ( 1 ) *0 .0 

BK ( 1 ) =0 • 0 

LM=1 

T=T I 

12 T=T*DT 

X*2.«EXP(-D*T)-1. 

CALL JACOBI (MN*XfB.P) 

SF<1)»CU>»PU> 

DO 10 L=2*MN 
L 1 =L“ 1 
AL=L 

SF(L)=SF(L1)*C(L)*°(D 
10 CONTINUE 
lm=lmm 

BK(LM)sSF(5) 

TT (LM)sT 

IF(T.LE.TF) 60 TO 12 
PRINT 97 

97 FORMAT (/////5Xf* T K T K 

IT K *) 

DO 31 MYsl.25 
MA=MY ♦ 1 
MB=MA>25 
mC=mB»25 
mD = mC ♦ 2*^ 

PRINT 96fTT(MA),BKrMA) fTT(MB) fBK(MB) ,TT(MC)fRK (MC) *T1 
96 FORMAT (5X f F5.2 ♦ 3X . ^7 . 5 * 3X f F5 • 2 « 3X f F7 .5 . 3X ,F5 .? , lx t F7 , 
1F7.S) 

31 CONTINUE 
RETURN 
END 


SYMBOLIC REFEPENCE MAP (P=l) 


vartarles 

SN TYPE 

1 

151 

AL 

PEAL 


4 

BK 

REAL 

ARRAY 

0 

D 

REAL 


147 

L 

INTEGER 


150 

LI 

INTEGER 


154 

MB 

INTEGER 


156 

MO 

INTEGER 


152 

MY 

INTEGER 


157 

SF 

REAL 

ARRAY 

1 

TF 

REAL 


151 

TT 

REAL 

ARRAY 


RELOCATION 


F.P. 


0 

B 

real 

0 

C 

PEAL 

2 

DT 

real 

1 44 

LM 

INTEGEP 

153 

ma 

integer 

155 

MC 

intfger 

3 

MN 

integer 

241 

P 

real 

145 

T 

real 

0 

T I 

peal 

146 

X 

peal 


1 


5 


10 


15 


20 


25 


SUBROUTINE JACOB I <N * X * B.PB ) 

C THIS PROGRAM CALCULATES JACOBI POLYNOMIALS OF OROER 

C K-l WITH ARG X AND PARAMETER B GT -1 

DIMENSION PB(N) 
an=n 

IF ( AN-2 • ) 1 « 2 ♦ 3 

1 PB ( 1 ) *1 • 

RETURN 

2 PB(l)sl. 

PB(2)=X-B*(l.-X)/2. 

RETURN 

3 BSOsB*B 
BONE*B* 1 • 

PB<1)=1. 

PB(2)*X-B*(l.-X)/2. 

00 4 K = 3 ♦ N 

ak=k 

AK 1 =AK-1 . 

AK2=AK-2. 

Kl=K-l 
K2=K-2 

C01*( (2.*AK1)*R)*X 
C01*( (2.*AK2)*B)*Cnl 
C01 = ( <2.*AK2> ♦90NE1 « <C01-9SQ> 

C02s2.*AK2« (AK2+8) * ( (2.*AK1 ) *8) 

C0s2.*AKl»(AKl+R)*f (?.*AK2> *9) \ 

4 PB(K)=(C01*PB(Kl)-r0?»PB(K2) )/C0 i 

RETURN ] 

END "i 

•j 

j 

( 


SYMBOLIC REFERENCE MAP <P=!) 

ENTRY POINTS 
3 JACOB! 


VARIABLES SN 

TYPE 

RELOCATION 





10* 

AK 

REAL 




107 

AK1 

real 

1 1 0 

AK? 

PEAL 




102 

AN 

PEAL 

0 

8 

REAL 


F.P, 


104 

BONE 

peal 

103 

BSO 

PEAL 




115 

CO 

PEAL 

113 

C01 

REAL 




114 

CC2 

REAL 

105 

K 

INTEGER 




111 

K 1 

INTEGER 

112 

K2 

INTEGER 




0 

N 

INTEGER 

0 

PB 

REAL 

ARRAY 

F.P. 


0 

X 

real 

STATEMENT LABELS 








0 

1 

INACTIVE 


24 

2 




0 

u 








LOOPS 

LABEL 

INDEX FROM-TO 

length 

PROPERTIES 


47 

4 

K 

16 27 

2RB 


OPT 




i 


1 

5 


10 


15 


20 


25 


30 


35 


40 


45 


SYMBOL i ' 


SUBROUTINE SPLINE (y*Y*M*CfXINT*YINT) 

DIMENSION X<50> *Y(S0) *C<4*50) 

IF (XINT-X ( 1 ) ) 1 * 1 0 * 1 1 

10 Y INT*Y ( 1 ) 

RETURN 

11 CONTINUE 
IF(X(M)-XINT)1» 12*13 

12 YINT*Y(M) 

RETURN 

13 CONTINUE 
K*M/2 
N=M 

2 CONTINUE 

IF (X(K)-XINT)3»14» C 

14 YlNTsY(K) 

RETURN 

3 CONTINUE 

IF < XINT“X (K ♦ 1 ) ) 4 * 1 c * 7 

15 Y INT=Y ( K ♦ 1 ) 

RETURN 

4 CONTINUE 

YINT=(X<K*1 J-XINT)*(C(1 « K > * < X (K ♦ 1 > -X INT > «*2*C < 3* K ) ) 
YINT=Y1NT* (XINT-X {*) )*(C(2*K)*(XINT-X(K) )#*2*C(4*K) ) 
RETURN 

5 CONTINUE 

IF(X (K-1)-XINT>6»14,17 

6 K*K-1 
GO TO 4 

16 YINT=Y(K-1) 

RETURN 

17 N=K 
K=K/2 
GO TO 2 

7 LL=K 

K* (N*K) /2 

8 CONTINUE 

IF(X(K)-XINT)3*14,)8 

18 CONTINUE 

1F(X(K-1)-XINT)6* 16*19 

19 N=K 
K=(LL*K)/2 
GO TO 8 

1 PRINT 101 

101 FORMAT <* OUT OF RANGE FOR INTERPOLATION *) 

STOP 

ENO 


REFERENCE MAP (R=l) 


ENTRY POINTS 
3 SPLINE 


1 SUBROUTINE SPL ICE < X * Y *M,C ) 

DIMENSION X<50) *Y<50) *0(50) *P(50) *E(S0> *C(4*50> 
DIMENSION A (50*3) ♦ R (50 ) • Z (50 ) 

MMsM-1 

5 00 2 K*1*MM 

0(K)«X(K4l)-X(K) 

P(K)*0(K)/6. 

2 E<K)*(Y(K*1)-Y(K) ) /D(K> 

DO 3 K*2*MM 

10 3 0<K)aE(K)-E<K-l) 

A(l*2)=-1.-D(l)/D(?) 

A ( 1 *3) »D ( 1 ) /D ( 2) 

A<2*3)=P(2)-P(1)*A<1,3) 

Al2,2)=2.»IPIlhP|?M-P(l)*A|1.2) 

15 A(2*3)sA(2*3)/A(2«?) 

0(2) =8 (2) /A (2*2) 

DO 4 K=3*MM 

A <K*2)*2.*(P(K-1 ) ♦o(K) )-P(K-l > »A (K-l *3) 
0(K)=B(K)-P(K-1 ) *0 fK-1 ) 

20 A(K,3)=P(K)/A(K,2) 

4 B(K)s0(K)/A(K*2) 

0=D(M-2)/0(M-l ) 

A(M*i)*l.»0*A(M-2,D 
A (M*2 ) s-O-A (M*1)*A(M-1*3) 

25 B(M)-0(M-2)-A(M,l)«B(M-l) 

Z(M)=B(M)/A(M*2) 

MN=M-2 

DO 6 I=1*MN 

K=M-I 

30 6 Z(K)=B(K)-A(K*3)»7fK*l> 

Z(1)*-A(1»2)«Z(?)-A (1 *3)*ZC3) 

DO 7 K=1 *MM 
Q=1*/(6.*D<K) ) 

C(1.K)=Z(K)*Q 

35 C(2*K)=Z(K*1)*Q 

C(3*K)=Y (K)/D(K)-Z(K)*P(K) 

7 C(4*K)=Y ( K ♦ 1 )/D(K)-Z (K*1)*P(K) 

RETURN 

END 


SYMBOLIC REFERENCE MAP (R=l > 

F.NTPY POINTS 
3 SPLICE 


VARIABLES 5N TYPE RELOCATION 


373 

A 

REAL 

ARRAY 


621 

R 

REAL 

0 

C 

REAL 

ARRAY 

F »P . 

145 

0 

real 

31 1 

E 

REAL 

ARRAY 


144 

1 

INTEGER 

141 

K 

INTEGER 



0 

M 

INTEGER 

140 

MM 

INTEGER 



143 

MN 

INTEGER 

227 

P 

REAL 

ARRAY 


142 

Q 

REAL 

0 

X 

PEAL 

ARRAY 

F • P . 

0 

Y 

REAL 


